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It has been known that epidemic outbreaks in the SIR model on networks are described by 
phase transitions. Despite the similarity with percolation transitions, whether an epidemic outbreak 
occurs or not cannot be predicted with probability one in the thermodynamic limit. We elucidate 
its mechanism by deriving a simple Langevin equation that captures an essential aspect of the 
phenomenon. We also calculate the probability of epidemic outbreaks near the transition point. 
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I. INTRODUCTION 

We start with the following question: How can it be 
determined whether an epidemic outbreak has occurred. 
Obviously, this is hard to answer, because an accurate 
model of epidemic spread in real societies, which in- 
clude complicated and heterogeneous human-to-human 
contact, cannot be constructed. Then, is it possible to 
predict the outbreak for a simple mathematical model? 
Even in this case, the manner of the early spread of dis- 
ease may significantly influence states that manifest after 
a sufficiently long time. For example, it seems reasonable 
to conjecture that whether a single infected individual 
with a very high infection rate causes an outbreak may 
depend on the number of people infected by the individ- 
ual, which is essentially stochastic. In the present paper, 
we attempt to formulate this conjecture. 

Specifically, we study the stochastic SIR model as the 
simplest epidemic model, where an edge in the network 
represents a human-to-human contact and the infection 
rate A (the infection probability per unit time in each 
edge) is a parameter of the SIR model (see e.g. Ref. [l| 
for an introduction to the stochastic SIR model; see also 
Refs. 0, |3j for related social dynamics on complex net- 
works). The SIR model may be defined for well- mixed 
cases homogeneous networks and scale- 

free networks |15l - ll8l | . A remarkable phenomenon is that 
when A exceeds a critical value Ac, a disease spreads 
to macroscopic scales from a single infected individual, 
which corresponds to an epidemic outbreak. This was 
found in well-mixed cases and random graphs, but Ac = 
for scale free networks. That is, epidemic outbreaks are 
described as phase transition phenomena. In addition 
to the interest in theoretical problems, recently, the SIR 
model on networks has been studied so as to identify in- 
fluential spreaders |19|] and so as to determine a better 
immunization strategy [l^, [2l| . 

Although the phase transition in the SIR model may 
be a sort of percolation transition, its property is differ- 
ent from that of standard percolation models. In the SIR 
model exhibiting the phase transition, the order param- 
eter characterizing it may be the fraction of the infected 
population, which is denoted by p. Indeed, p = in the 
non-outbreak phase (A < Ac), whereas the expectation 



of p becomes continuously non-zero from when A > Ac. 
This phenomenon is in accordance with the standard per- 
colation transition. However, on one hand, the order pa- 
rameter in the percolated phase, e.g. the fraction of the 
largest cluster, takes a definite value with probability one 
in the thermodynamic limit; on the other hand, the frac- 
tion of the infected population in the SIR model is not 
uniquely determined even in the thermodynamic limit. 
In fact, it has been reported that the distribution func- 
tion of the order parameter in SIR models with finite sizes 
shows two peaks at p = and p = p* for well-mixed cases 
, hornogeneous networks IHI [3 1 ^^^^ scale- free 
networks [18'| . Mathematically, the probability density of 
p in the thermodynamic limit may be expressed as 

P{p- A) = (1 - q{\mp) + qimp - p.), (1) 

where q — for A < Ac and q ^ for A > Ac . This means 
that the value of the fraction of the infected population 
in the outbreak phase, which is either or p*(A), cannot 
be predicted with certainty. We call this phenomenon 
the intrinsic unpredictability of epidemic outbreaks. 

In this paper, we clearify the meaning of (P). We first 
observe the phenomenon in the SIR model defined on a 
random regular graph. By employing a mean field ap- 
proximation, we describe the epidemic spread dynamics 
in terms of a master equation for two variables. Then, 
with a system size expansion, we approximate the so- 
lutions to the master equation by those to a Langevin 
equation. Now we can analyze this Langevin equation 
and work out the mechanism of the appearance of the 
two peaks. We also calculate q{\) near the transition 
point. 



II. MODEL 

Let G be a random fc-regular graph consisting of N 
nodes. For each x £ G, the state a{x) S {S,I, R} is 
defined, where S, I, and R represent Susceptible, Infec- 
tive, and Recovered, respectively. The state of the whole 
system is given by {ax)xGG: which is denoted by cr col- 
lectively. The SIR model on networks is described by 
a continuous time Markov process with infection rate 
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A and recovery rate Concretely, the transition rate 
T4^(cr —> cr') of the Markov process is given as 



with 



Wicr -^cr') = Y^ w{(T ct'\x), 



w{cr — )■ (t'\x) 



(2) 



yeB{x) 



(3) 




where B{x) is a set of fc-adjacent nodes to x £ G. Here- 
inafter, without loss of generality, we use dimensionless 
time by setting fi — 1. For almost all time sequences, 
infective nodes vanish after a sufficiently long time, and 
then the system reaches a stationary state, which is called 
the final state. The ratio of the total number of recovered 
nodes to N in the final state is equivalent to the fraction 
of the infected population p. This quantity measures the 
extent of the epidemic spread. At t = 0, we assume that 
a = I for only one node selected randomly and that a = S 
for the other nodes. 

In Fig. [Tl as an example, we show the result of numer- 
ical simulations for the model with fc = 3 and N = 8192. 
We measured the probability density P{p; A) of the frac- 
tion of the infected population p for various values of A. 
This figure suggests that the expectation of p becomes 
non-zero when A exceeds a critical value. The impor- 
tant observation here is that log P in the outbreak phase 
has a sharp peak near p — 0, too. Indeed, the inset 
in Fig. [T] clearly shows the existence of the two peaks 
in logP with A = 1.5. Similar graphs were reported 
in Refs. [1, 0, 1 [H [ll IHl • The existence of the two 
peaks is not due to a finite size effect, as shown in Fig. [21 
where the probability that p > 1/16, which is denoted by 
p{p > 1/16), is plotted as a function of A for several val- 
ues of N. Note that limjv->ooP(p > 1/16) = q{X) when 
p*{X) > 1/16. These results suggest the limiting density 
([T]), where q{X) becomes continuously non-zero for A > Ac 
whereas q{X) = for A < Ac. This is the phenomenon 
that we attempt to understand in this paper. 



III. ANALYSIS 

Defining two variables s = '^^Slax, S)/N and i = 
X^a; ^(''^^' -^)/-^' consider a continuous-time Markov 
process of the two variables as an app roximation of 
the SIR model on the network [2^, [23]. We expect 
the phenomenon we are concerned with to be repro- 
duced within this approximation; we verify this at a 
later stage. The transition rate of (s,«) — >■ {s,i — 1/N) 
is exactly given as Ni, and we approximate the rate 
(s,«) — ^ (s — l/N,i + 1/N) as XkNsip, where ip is the 
probability of finding y e B{x) such that Cy = I for any 
X. Here, the infective nodes form a connected cluster, 
and this cluster is tree-like because the typical size of the 



FIG. 1: (Color online) Color representation of logP(p;A) in 
the (p. A) plane. It is obtained by numerical simulations of 
the SIR model on a random regular graph. The inset shows 
logP(p; A) as a function of p for A = 1.5. 
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FIG. 2: (Color online) p{p > 1/16) as a function of A for 
several values of A'^. 



loops is 0(log N). Now, as an approximation, we assume 
that there are Ni{k — 2) edges connecting the tree-like 
cluster with susceptible nodes |2^,|25|. Therefore, ip is es- 
timated as the rate of Ni{k—2) to the number of all edges 
Nk in the thermodynamic limit. That is, -0 = i{k — 2)/fc. 
Below, we focus on the case fc = 3. 

Let P{s, i, t) be the probability density of s{t) = s and 
i{t) — i. Then, P{s^i,t) obeys the master equation 



dP{s,i,t) 
dt 



= N[i + 



+ NX[s + 



1 

iV 

NXsiP {s,i,t) 



1 

iV 



P 



1 

s H ,i 



1 



(4) 



When N is sufficiently large, the master equation for 
P(s, i, t) can be expanded as 



^+d,J, + dsJs + o(j^]=0, (5) 
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FIG. 3: (Color online) Result of numerical simulations of 
(O and (|8]). The presentation is the same as those in Fig. (2] 
= 1024 for the main frame and = 8192 for the inset. 



FIG. 4: (Color online) Functional forms of U{Y). X = 0.5 
(red solid line) and A = 1.2 (green dotted line.) 



with 



(As - l)iP 
Js = —XsiP — ds 



Xsi 
2N' 



(As + l)i 
2N 

P]+^^ 



P 

Xsi 



-P 



(6) 



By assuming that 0{1/N'^) terms can be ignored, we 
obtain the Fokker-Planck equation [26] . 

It can be confirmed by direct calculation that this 
Fokker-Planck equation ([5]) describes the time evolution 
of the probability density for the following set of Langevin 
equations: 



ds 
'dt 
di 
di 



'Xsi 



Xsi 
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(7) 

(8) 



where is Gaussian white noise that satisfies {^i (t)) =0 
and {^i (t) {t')) = SijS {t - t'). The symbol "•" in front 
of ^1 and ^2 in (O and ^ represents the Ito product 
rule. The same equations as ([7]) and dS]) were presented 
in Refs. [2^[2^. In this description, the fraction of the 
infected population is given by 



p — 1 — s(cxd). 



(9) 



In Fig. [31 we show the result of numerical simulations of 
the Langevin equations ([7]) and ([8]). Comparing Fig. [3] 
with Fig. [U we find that the phenomenon under study is 
described by the Langevin equations ([7]) and ([5]). Thus, 
our problem may be solved by analyzing them. 

Now, the key idea of our analysis is the introduction 
of a new variable Y = y/iN. Then, ([7|) and dH) are re- 
written as 



-AsF^ _ VAsy2 . 



ds _ 1 

~ TV 

^ ^ i|(As-l)y-i(As + l)l 

dt 2 \^ ' 4,^ ' Y 

1 



(10) 



(11) 



where it should be noted that the multiplication of the 
variable Y and the noise does not appear in ([TT|). We 
then consider the probability q{X) in the thermodynamic 
limit as the probability of observing Y ~ N^/'^, because 
it is equivalent to p > 0. 

Here, from (fTU)) and (fTTj) . we find that the character- 
istic time scale of s is iV times that of Y . Thus, when 
iV is sufficiently large, s almost retains its value when Y 
changes over time. In particular, it is reasonable to set 
s = 1 when t is shorter than N . In this time interval, 
PT|) is expressed as 



where D 
as 



'2D^, (12) 
(A+l)/8 and the potential U{Y) is calculated 



U{Y) = -1 (A - 1) + i (A + 1) log(y). (13) 

f is Gaussian white noise with unit variance, where we 
have used the relation VA/2^i-h 1/2^2 = \/A + 1/2^. The 
initial condition is given as F(0) = 1. It should be noted 
that ((T^ is independent of N . Thus, solutions satisfy- 
ing Y ~ in (|10p and (|TT|) correspond to solutions 
satisfying y — > cx) in ()12|) . We identify q{X) with the 
probability of finding these solutions. We now derive 
this probability. 

First, we investigate the shapes of the graph U {Y). 
We find that C/(0+) = -oo for any A and that U {Y) 
monotonically increases in F for A < 1, while U{Y) has 
a single maximum peak at 1" = for A > 1, where 



K = 



1 /A-1- 1 



A- r 



(14) 



As a reference, in Fig. 21 we show the shapes of U {Y) 
for A = 0.5 and 1.2. 

Next, based on the shapes of the potential function, 
we discuss the expected behavior of solutions to (|12p . 
When A < 1, the probability of 1^ — cx) is obviously zero 
because U{Y) is a monotonically increasing function in 
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Y. That is, q{X) = in this case. The behavior for 
A > 1 is comphcated. We thus focus on the case that 
A = 1 + e, where e is a small positive number. In this 
case, ~ e~^/^. We then note that if a solution Y to 
p2)) happens to exceed , it is comparatively likely that 
y — > oo. Assuming that the probability of y — > oo under 
the condition y > y* at some time is unity, we estimate 
q{X) as the probability that Y exceeds y*. Furthermore, 
we express q{X) in terms of the transition rate T from 
y = 1 to y = y*. Noting that the transition rate from 
y = 1 to y = is equal to the recovery rate in the 
original SIR model, we can write 



T 



q = 



1 + T 



(15) 



Since T is positive and finite, we obtain < q{X) < 1. In 
this manner, we have clearly explained the probabilistic 
nature in the outbreak phase, and we have obtained Ac = 
1. 

Finally, we calculate q{X) quantitatively near the tran- 
sition point. From y, ~ e~^/^ and ?7(y,) ~ loge, we 
estimate the slope of the straight line connecting two 
points (l,t/(l)) and (y»,C/(y»)) in the {Y,U) plane as 
(t/(y*) - t/(l))/(y, - 1) ~ Ve(loge), which approaches 
zero in the limit e 0. Thus, the transition from y = 1 
to y = y* may be assumed to be free Brownian mo- 
tion with the diffusion constant D — {X + l)/8. The 
transition rate from y = 1 to y* is then estimated as 
T = 2D/Y^^ = € + 0{e^). We thus obtain 



g(A) = e + Oie'). 



(16) 



In Fig. [5l we compare the theoretical result with those 
obtained in numerical simulations of ((TU)) and dTTI) . We 
measured the probability that p > 0.003, which is de- 
noted as p{p > 0.003). Recall that liuiN^oo p{p > 
0.003) = q{X) when p^X) > 0.003. Since the experi- 
mental result suggests p{p > 0.003) = e + O(e^) in the 
limit N ^ oo, we claim that the theoretical result ([16]) 
is in good agreement with the experimental result. 



IV. CONCLUDING REMARKS 

In this paper, we have achieved a novel understanding 
of the intrinsic unpredictability of epidemic outbreaks by 
analyzing the Langevin equation (1121) . which effectively 
describes this singular phenomenon. Further, trajecto- 
ries in the outbreak phase are divided into two groups: 
trajectories in one group are absorbed into zero, and the 
others diverge in (fT^ . The division corresponds to the 
non-trivial limiting density given in ([T]). On the basis of 
this description, we calculated the probability of an epi- 
demic outbreak near the transition point. Before ending 
the paper, we make a few remarks. 

First, the probability o( A) was studied in the mathe- 
matical literature (see [l^ and [l^ as reviews.) To the 
best of our knowledge, the method proposed in this pa- 
per has never been used in previous studies. It might 
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FIG. 5: (Color online) p{p > 0.003) as a function of e = A — 1 
obtained by numerical simulations of (|10p and The guide 
line represents p{p > 0.003) = e, which is expected from the 
theoretical analysis. 
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FIG. 6: Probability that p > 1/16 in the SIR model on a 
random regular graph; m — N/128. 



be interesting to connect our analysis with mathematical 
studies. 

Second, although we have investigated the simplest 
model in this paper, similar analysis might be applied 
to various models. For example, we can consider the 
case that there are m infected nodes at time t = 0. Since 
the essence of the phenomenon is the existence of y*, 
the same result is obtained when m is independent of 
N. However, for the case m = cN with a small positive 
number c, Y{t) is never adsorbed to zero in the outbreak 
phase, because y(0) is infinitely far away from Y = Y^,. 
This is qualitatively different from the case m — 1, which 
was reported in Refs. [29l, fsol]. In fact, as suggested in 
Fig. [6l q{X) jumps discontinuously to q{X) = 1 which is 
similar to the behavior observed in standard percolation 
transitions. 

Finally, as another generalization, one may study the 
behavior of the SIR on more complex networks. In these 
cases, since the mean field approximation might not be 
effective, one needs to devise a new technique to describe 
the unpredictability of outbreaks. Moreover, one of the 
most interesting is to predict probabilistic epidemic out- 
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breaks from limited data on realistic networks. We hope 
that future studies will address these problems. 
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